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Abstract 



We review a recently devised Monte Carlo simulation method for the 
direct study of quasi-stationary properties of stochastic processes with 
an absorbing state. The method is used to determine the static cor- 
relation function and the interparticle gap-length distribution in the 
critical one-dimensional contact process. We also find evidence for 
power-law decay of the interparticle distance distribution in the two- 
particle subspace. 
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I. INTRODUCTION 



Stochastic processes with an absorbing state arise frequently in statistical physics 
[1,2], epidemiology [3] and related fields. Phase transitions to an absorbing state in 
spatially extended systems, exemplified by the contact process [4,5], are currently 
of great interest in connection with self-organized criticality [6], the transition to 
turbulence [7], and issues of universality in nonequilibrium critical phenomena [8-10]. 

The quasi- stationary (QS) distribution, (that is, conditioned on survival), is 
very useful in the study of processes with an absorbing state model. Conventional 
simulations of "stationary" properties of lattice models with an absorbing state 
actually study the quasi-stationary regime, given that the only true stationary state 
for a finite system is the absorbing one. We recently devised a simulation method 
that yields quasi-stationary properties directly [11]; it samples the QS probability 
distribution just as conventional Monte Carlo simulation samples the equilibrium 
distribution. Here we use the method to study the static correlation function and 
other configurational properties of the critical contact process, the prime example 
of an absorbing-state phase transition. 

In the following section we review the basis of our method. Then in Sec. Ill we 
apply it to determine the static two-point correlation function of the contact process 
on a ring. We summarize our findings in Sec. IV. 



II. BACKGROUND 



Consider a continuous-time Markov process X t taking values n = 0, 1, 2, S : 
with state n = absorbing. We use p n (t) to denote the probability that X t = n, given 
some initial state X . The survival probability is P s {t) = J2n>iPn{t) = 1 — Po{t)- 
We suppose that as t — > oo the p n , normalized by the survival probability P s (t), 
attain a time-independent form, thus defining the quasi-stationary distribution p n : 

with p = 0. (We further assume that the limiting distribution does not depend on 
the initial state, as long as X ^ 0.)The QS distribution is normalized so: 

EPn=l- (2) 

n>l 



As shown in [12], the QS distribution is the stationary solution to the following 
equation of motion (for n > 0) 
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dq n 

— = -w n q n + r n + r q n , 



(3) 



where w n = J2 m w m,n is the total rate of transitions out of state n, and r n = 
Hm w n,mQm is the flux of probability into this state. To see this, consider the master 
equation (Eq. (3) without the final term) in the QS regime. Substituting q n (t) = 
P s (t)p n , and noting that in the QS regime dP s /dt = — r = — P s J2m w o,mP m , we see 
that the r.h.s. of Eq. (3) is identically zero if q n = p n for n > 1. The final term in Eq. 
(3) represents a redistribution of the probability tq (transfered to the absorbing state 
in the original master equation), to the nonabsorbing subspace. Each nonabsorbing 
state receives a share equal to its QS probability. 

In [11] we introduce a process X*, whose stationary probability distribution is 
the quasi- stationary distribution of X t . (Note that in order to have a nontrivial 
stationary distribution, X 4 * cannot possess an absorbing state.) The probability 
distribution of X* is governed by Eq. (3), which implies that for n > (i.e., away 
from the absorbing state), the evolution of X t * is identical to that of X t . When X t 
enters the absorbing state, however, X t * instead jumps to a nonabsorbing one, and 
then resumes its "usual" evolution (with the same transition probabilities as X t ), 
until such time as another visit to the absorbing state is imminent. 

In Eq. (3) the distribution q n is used to determine the value of X 4 * when X t visits 
the absorbing state. Although one has no prior knowledge of q n (or its long-time 
limit, the QS distribution p n ), one can, in a simulation, use the history X* (0 < s < 
t) up to time t, to estimate the q n . This is done by saving a sample {n\,n 2 , ...,um} 
of configurations visited. We update the sample by replacing from time to time 
one of the configurations with the current one. In this way the distribution for the 
process X t * will converge to the QS distribution (i.e., the stationary solution of Eq. 
(3)) at long times. Summarizing, the process X t * has the same dynamics as X u 
except when a transition to the absorbing state is imminent: XI then is placed 
in a nonabsorbing state, selected at random from a sample over the history of the 
realization. (In the simulation, a list of M configurations is maintained. Whenever 
the time increases by 1, the list is updated with probability p rep , by replacing a 
randomly chosen configuration on the list with the current one.) In effect, the final 
term in Eq. (3) is represented as a memory in the simulation. 

The above scheme was shown [11] to yield precise results, in accord with the 
exact QS distribution for the contact process on a complete graph [12], and with 
conventional simulations of the same model on a ring, for which exact results are not 
available. QS simulation results were found rather insensitive to the choice of list 
size M and replacement rate p rep . In the studies reported below we use M = 1000 
and p re p = 0.001. 
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III. CONTACT PROCESS: CORRELATION FUNCTION 



The contact process (CP) [4] is a continuous-time Markov process on a lattice, 
in which each site % is either occupied = 1), or vacant (cx;(t) = 0). Transitions 

from <7j = 1 to Oi — occur at a rate of unity, independent of the neighboring 
sites. The reverse transition can only occur if at least one neighbor is occupied: the 
transition from <j{ — to Oi — 1 occurs at rate Ar, where r is the fraction of nearest 
neighbors of site i that are occupied; thus the state Oi = for all i is absorbing. (A is 
a control parameter governing the rate of spread of activity.) The order parameter 
p = (o"j) is the fraction of occupied sites. The model exhibits a continuous phase 
transition at A c = 3.297848(20) [13]. For A < A c , the stationary value of p is zero. 

The CP has attracted much interest as a prototype of a nonequilibrium critical 
point, a simple representative of the directed percolation (DP) universality class. 
Since its scaling properties have been discussed extensively [8-10] we review them 
only briefly. As the critical point is approached, the correlation length £ and cor- 
relation time r diverge, following £ oc |A| _1/± and r oc |A| -I/ ii, where A = A — A c 
is the distance from the critical point. The order parameter scales as p oc A 13 for 
A > 0. At the critical point the quasi-stationary value of the order parameter scales 
as: p oc L-M v± . 

An aspect of the CP that has not, to our knowledge, been studied in simulations 
is the static correlation function. Of particular interest is how correlations decay at 
the critical point. The correlation function is defined via 



where the average is over the stationary distribution of the process. 

Based on experience with equilibrium critical phenomena, we expect C(r) to 
decay as a power law at the critical point. The decay exponent can be determined 
via a scaling argument. To begin, we normalize C(r) to its value at r = 0: 



(L denotes the lattice size.) At the critical point, x ~ L 1 ^ 1 - with 7 = du± — 2j3 [8]. 
A simple calculation yields 




(4) 




(5) 



Now consider the scaled variance 




(6) 




(7) 
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where in the last step we used translation invariance, and assumed a hypercubic 
lattice of L d sites. Now suppose c(r) ~ r~ a for large r. Approximating the sum by 
an integral, and recalling that p ~ L^ 13 ^ 1 - at the critical point, we find 

X ~ L^p/u^ _ L d-a-/3/u ± ^ 

implying that a = j3/vj_ and 

C(r) ~ pc(r) oc (rL)-^ v± (9) 

in the critical stationary state. The relation a = j3/vj_ was demonstrated some time 
ago by Grassberger and de la Torre, who showed that C(r) ~ r - 25 / z a t the critical 
point [14]. (Their result is seen to be equivalent to ours when we recall the scaling 
relations z = 2v±/v\\ and 5 = f3/v\\-) 

We use the QS simulation method outlined above to determine the correlation 
function. The process is simulated in five independent realizations of 2 x 10 8 time 
steps. As is usual, annihilation events are chosen with probability 1/(1 + A) and 
creation with probability A/(l + A). A site is chosen from a list of currently occupied 
sites, and, in the case of annihilation, is vacated, while, for creation events, a nearest- 
neighbor site is selected at random and, if it is currently vacant, it becomes occupied. 
The time increment associated with each event is At = l/N occ , where N occ is the 
number of occupied sites just prior to the attempted transition [8] . 

In Fig. 1 we plot C*(r) = L^C(r) for L = 1280 and 2560, using the best 
available estimate (from series analysis), /3/v± = 0.252072(8) [15]. The data collapse 
for the two lattice sizes is nearly perfect. For r L the correlation function 
indeed follows a power law C* ~ r -i 3 / u ±^ while for r = L/2 it attains a minimum, 
as expected due to the periodic boundaries. To determine the decay exponent 
we analyze the local slope a(r) (see Fig. 2), obtained from a linear fit to the 
data for InC* versus lnr, using points equally spaced in lnr, in finite intervals 
[r ,3r ]. For r <C L the local slope is nearly constant, but it of course veers 
upward as r approaches L/2. We therefore perform an extrapolation (to r — > oo) 
of the local slope versus 1/r, using only the data on which the results for the two 
lattice sizes agree, to eliminate finite-size effects. The result of this extrapolation is 
(3/v± = 0.252(1), consistent with the best estimate. 

For A < A c , the correlation function decays exponentially, C(r) ~ e _r ^. This 
is evident in the inset of Fig. 1, where we plot C = r^^ ± C(r). Exponential decay 
is clear for A = 0.99A C ; linear regression yields £ ~ 356 in this case, well below the 
system size (L = 2560) in this study. Interestingly, the decay of C is also perceptible 
for A = 0.999A C , though partly masked due to the finite system size. In general the 
decay of correlations should be evident for L > £, where £ is the correlation length 
in the infinite-size limit. Since £ ~ |A| _1/± , deviations from criticality on the order 
of |A| ~ L~ l l Vl - (or greater) should be detectable in the correlation function for a 
system of size L. 



5 



In the critical stationary state, the distribution of particles is scale-invariant, as 
reflected in the power-law decay of C(r). The distribution of gaps, or strings of 
empty sites between successive occupied sites also follows a power law. (A gap of 
size g corresponds to sites i and i + g + 1 occupied and all intervening sites empty.) 
We determine the gap-size distribution P(g) (normalized to the number of gaps in 
the configuration). As shown in Fig. 3, P(g) exhibits a power law decay, P ~ g~ T , 
with t ~ 1.70, over an intermediate range that appears to grow with system size. 
Analysis of the local slope yields r = 1.73(1). 

We can relate the exponent r to other critical exponents via a simple scaling 
argument. Since there is one gap per particle, the mean gap size (g) is just the 
reciprocal of the particle density. Thus in a system of size L at the critical point, 
(g) ~ L /3 ^ u± . Assuming P(g) ~ g~ T for g > 1, we have 



implying r = 2 — (3/v± ~ 1.748. Our simulation result is about 1% smaller than 
the value predicted by the scaling argument. The discrepancy is likely caused by 
finite-size corrections that limit the power-law regime of P(g)- 

Finally, we report preliminary results on a surprising behavior of the interparticle 
distance in the two-particle subspace. Let d denote the separation between the 
occupied sites, given that that exactly two sites are occupied. (We take the minimum 
distance under periodic boundaries, so that d < L/2.) Since particles are highly 
clustered in the critical CP, we should expect the two-particle distance distribution 
P2(d) to decay with separation. Our results (see Fig. 4) from QS simulations at 
the critical point suggest a power-law decay, P2(d) ~ d~ K , with k ~ 2.45. We have 
no way of relating this exponent to the known critical exponents. Indeed, whether 
P2{d) follows a power-law will have to be confirmed in larger-scale simulations. This 
is somewhat challenging since, as the system size increases, the probability of having 
exactly two particles becomes ever smaller. 



We have applied a new simulation method for quasi- stationary properties to 
determine the static correlation function C(r) of the critical contact process in one 
dimension. Our results support the behavior C(r) ~ l/(rL) /3 / iy± anticipated from 
scaling arguments. We also studied the gap-size distribution, which shows evidence 
of a power-law decay, P(g) ~ g~ T , with r ~ 2 — j3/u±, the value predicted by 
scaling. Finally we note an apparent scale-invariant behavior of the interparticle 
distance distribution P2^d) in the two-particle subspace. Study of the correlation 
function and the gap-size distribution promise to be useful in characterizing scaling 
behavior of new models, and may also be useful in locating the critical point. 




(10) 



IV. SUMMARY 
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FIGURE CAPTIONS 



FIG. I. QS simulation results for the scaled correlation function C* = L l3 ^ u± C(r) in 
the critical one-dimensional contact process. Symbols: x: L = 1280; +: L = 2560. 
The slope of the straight line is -0.252. Inset: semi-logarithmic plot of C = r l3 ^ u± C(r) 
versus r in a system of 2560 sites, for A = A c (upper curve), 0.1% below A c (middle) 
and 1% below A c (lower). 

FIG. 2. Local slope —ct(r) of the correlation function versus 1/r. Open symbols: 
L = 1280; filled symbols: L = 2560. 

FIG. 3. Gap-length distribution in the critical one-dimensional CP. □: L — 640; +: 
L = 5120. The slope of the straight line is -1.73. 

FIG. 4. Distribution of interparticle distances d in the two-particle subspace of the 
critical CP. Open symbols: L = 640; filled symbols: L = 1280. The slope of the 
straight line is -2.45. 
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